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Abstract 

In this paper we develop new techniques for revealing geometrical structures in 
phase space that are valid for aperiodically time dependent dynamical systems, which 
we refer to as Lagrangian descriptors. These quantities are based on the integration, 
for a finite time, along trajectories of an intrinsic bounded, positive geometrical and/or 
physical property of the dynamical system. We discuss a general methodology for con- 
structing Lagrangian descriptors, and we discuss a "heuristic argument" that explains 
why this method is successful for revealing geometrical structures in the phase space 
of a dynamical system. We support this argument by explicit calculations on a bench- 
mark problem having a hyperbolic fixed point with stable and unstable manifolds that 
are known analytically. Several other benchmark examples are considered that allow us 
the assess the performance of Lagrangian descriptors in revealing invariant tori and re- 
gions of shear. Throughout the paper "side-by-side" comparisons of the performance of 
Lagrangian descriptors with both finite time Lyapunov exponents (FTLEs) and finite 
time averages of certain components of the vector field ("time averages") are carried 
out and discussed. In all cases Lagrangian descriptors are shown to be both more 
accurate and computationally efficient than these methods. We also perform compu- 
tations for an explicitly three dimensional, aperiodically time-dependent vector field 
and an aperiodically time dependent vector field defined as a data set. Comparisons 
with FTLEs and time averages for these examples are also carried out, with similar 
conclusions as for the benchmark examples. 

An understanding of the interrelation of geometrical structures in phase space 
and individual trajectories is central to the global, geometrical point of view 
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of dynamical systems theory, and over the past 30 years this point of view 
has played an important, and evolving, role in our understanding of transport 
in fluids from the Lagrangian point of view. In this paper we further de- 
velop a new technique, that we refer to as Lagrangian descriptors for revealing 
the geometrical structures in the phase space of a dynamical system, which 
also applies when the dynamical system has an arbitrary time dependence. 
Numerous examples are considered, from benchmark problems (where exact 
quantities can be computed) to examples at the forefront of research, such 
as a three dimensional Hills spherical vortex subjected to aperiodically time 
dependence and a velocity field defined as a data set obtained from the simu- 
lation of the wind driven, quasigeostrophic equations in a rectangular domain. 
For all examples "side- by-side" comparisons of the performance of Lagrangian 
descriptors with both finite time Lyapunov exponents (FTLEs) and finite time 
averages of certain components of the vector field ("time averages") are carried 
out and discussed. In all cases Lagrangian descriptors are shown to be both 
more accurate and computationally efficient than these methods. 



1 Introduction 



Since the insights of Poincare (1890) geometrical structures in phase space have played a 
central role in characterizing the global behaviours of dynamical systems. This point of 
view gives insight into the evolution of qualitatively distinct classes of trajectories without 
having to explicitly integrate and compare the behaviour of trajectories. This point of 
view has also been shown to provide much insight into fluid transport and mixing after 
the recognition that the equations for incompressible fluid particle motion (in the absence 
of molecular diffusion) are formally equivalent to Hamilton's equations where, in the fluid 
mechanics context, the stream function plays the role of the Hamiltonian function and the 
physical space in which the fluid moves plays the role of the phase space of the corresponding 



Hamilton's equations Aref (1984). This framework for studying fluid transport and mixing 



is often referred to as the "dynamical systems approach to Lagrangian transport" since 
the focus is on understanding the "organising structures" in phase space for fluid particle 
trajectories. References describing this approach for general fluid mechanics are Ottino 



(1989); Wiggins and Ottino 



Jones and Winkler (2002) 



(2006) 



( 2004 ) , and references specific to geophysical fluid dynamics are 



Wiggins (2005); Mancho et al. (2006c); Samelson and Wiggins 



In this paper we further develop a technique for revealing phase space structure that 



was proposed and used in Mendoza and Mancho (2010); Mendoza et al. (2010); Mendoza 



and Mancho (2012); de la Camara et al. (2012). This method is based on the computation 



of arc length of particle trajectories and it is extended by considering the integration of a 
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bounded, positive quantity that is an intrinsic geometrical and/or physical property of the 
dynamical system along trajectories of the dynamical system, for a finite time. Techniques 



based on the integration of properties along trajectories have been used in the past. Rypina 



et al. (2011 ) have proposed methods based on the idea of complexity of isolated trajectories. 



Recently in the computer graphics community, Acar et al. (2012) have also used the time- 
normalized arc length of trajectories for human action recognition. Their work is inspired 
on the line integral convolution discussed by Shi et al. (2008); Cabral and Leedom (1993) 
that advect conserved physical fields along trajectories. The advection of conserved physical 
fields has been also used in the atmospheric science community since the early 90s. |Q'Neill| 



et al. (1994); Sutton et al. (1994) proposed the reverse domain filling method for visualizing 



transport structures in geophysical flows. 

The techniques explained in the current article have been shown to be effective for 
revealing phase space structures for aperiodically time dependent flows. In Section [2] we 
describe the general methodology for constructing Lagrangian descriptors, and we empha- 
sise the role of the norm in quantifying the integral of the chosen positive quantity along 
trajectories. Section 2.1 describes the performance of Lagrangian descriptors in hyperbolic 
regions, and Section 2.1.1 discusses the benchmark problem of a linear, autonomous vec- 
tor field in the plane having a hyperbolic saddle point at the origin. This is an excellent 
example, not only for discussing the issues associate with the performance of Lagrangian 
descriptors, but also for comparing their performance with finite time Lyapunov exponents 
(FTLEs) and time averages of particular components of the vector field along trajectories 
since the stable and unstable manifolds of the hyperbolic point, as well as all trajectories, 
are known exactly. In this example we are able to illustrate the effect of the integration 
time and make precise the idea that "singular contours of the Lagrangian descriptors cor- 
respond to invariant manifolds", both numerically, and analytically in Section 2.1.2 This 
example also illustrates why the L 7 norms, 7 > 1 are not successful in revealing "sin- 
gular contours" in the same way as the L 1 norms, 7 < 1. In Section 2.3 we show that 
for this example FTLEs completely fail to reveal any phase space. Since the system is 
linear almost all FTLEs, for any time, are equal for every trajectory. A similar failure 
for finite time averages of a component of the vector field along trajectories is reported 
in Section 2.3. 2[ In Section |2.2| we consider another benchmark example-a linear vector 
field having an elliptic equilibrium point at the origin. As in the example of the linear 
saddle point, since this system is linear almost all FTLEs are zero, for any time. Hence 
the FTLE field fails to reveal any phase space structure. The Lagrangian descriptors, on 
the other hand, reveal the expected phase space structure, consistent with the trajectories. 
Since this example is linear elliptic point, effects of shear are not considered. In Section 



2.2.2 we consider an integrable example in action- angle coordinates (hence all FTLEs are 



zero, for any time) that not only contain regions of strong shear, but also invariant tori 
where the "twist condition" breaks down (i.e. "twistless tori"). We show that Lagrangian 
descriptors provide signatures of each of these features. An advantage of the linear saddle 
point and linear elliptic point examples is that they allow us to "isolate" hyperbolic and 
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elliptic behaviour. In Section [2.4| we consider the forced Duffing equation for three different 
types of time dependency: 1) no time dependency (forcing), i.e. the integrable case, 2) 
periodic time dependence, and 3) aperiodic time dependence (of a form that we define). In 
this example hyperbolic and elliptic behaviour are " intermingled" in a complicated manner 
(that we describe), and this affords us an ideal setting to compare Lagrangian descriptors, 
FTLEs, and averages of velocity components along trajectories " side- by-side" in different 
regions of the phase space of the forced Duffing oscillator. In Section [3] we consider a 
three dimensional vector field from fluid mechanics-the Hill's spherical vortex subjected 
to a time-aperiodic perturbation-that was studied in Branicki and Wiggins (2009). For 
this example we show that Lagrangian descriptors provide an efficient way of discovering 
and visualising detailed geometrical structures in the flow. A comparison is made with 
FTLEs and averages of components of the vector field along trajectories, and it is shown 
that these techniques may introduce "artefacts" that obscure the true geometric structures. 
In Section [4] we consider velocity fields defined as data sets. In particular, we consider a 
wind driven three-layer quasigeostrophic simulation in a rectangular domain in a "double 
gyre" configuration exhibiting aperiodic time dependence. This velocity field has been 
the subject of earlier studies (Coulliette and Wiggins (2001); Mancho et al. (2004, 2006c)) 
and therefore allows us to directly compare the performance of Lagrangian descriptors, 
FTLEs, and time averages of components of velocity along trajectories with distinguished 
hyperbolic trajectories (DHTs) and their stable and unstable manifolds. Additionally, we 
compare the convergence time of Lagrangian descriptors with FTLEs to the stable and 
unstable subspaces of a given DHT, where we see that the Lagrangian descriptors con- 
verge to these structures more quickly. In Section [5] we compare the computational effort 
required for Lagrangian descriptor computations with FTLE computations, and in Section 
[6] we summarise our conclusions. In three appendices we provide some additional technical 
details. 



Lagrangian Descriptors: Definitions, Heuristics, and An- 
alytical Results 



The concept of a Lagrangian descriptor was first introduced in Mendoza and Mancho 



(2010). There it was shown that this tool is able to provide a global dynamical picture 
of the geometric structures for arbitrary time dependent flows. In particular, Lagrangian 
descriptors are able to detect the main "organising centres" in the flow- hyperbolic tra- 
jectories and their stable and unstable manifolds and elliptic regions. In this section we 
describe precisely what we mean by the term "Lagrangian descriptor". We then provide 
a heuristic justification for why they "work". We build on this by then considering sev- 
eral benchmark examples where the phase space structure and stability characteristics are 
known analytically. This allows a deeper investigation of the heuristic justification, as well 
as a comparison with other commonly used method for extracting phase space structures, 
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such as finite time Lyapunov exponents (FTLEs). 

We consider a general time-dependent vector field on M. n : 

^ = v(x,t), xef.tel (l) 

where we assume that v(x, t) is C r (r > 1) in x and continuous in t. This is a sufficient 
condition for the existence of unique solutions that also allows for linearization. The time 
interval on which a solution exists is also an issue. However, we will proceed by assuming 
that the solutions exist for a sufficient time so that the integral expressions we derive have 
validity, and then we will verify the validity in specific examples. 



We begin by describing the Lagrangian descriptor first described in Mendoza and Man 



cho (2010). In that reference it was denoted by M, but henceforth we will denote it by Mi 
since we will introduce additional Lagrangian descriptors. M\ is the Euclidean arc length 
of the of the curve in phase space defined by a trajectory of that starts at x* at time 
t = t* for the time interval [t* — t, t* +r], i.e. 



Mi(k-,Ov,r = r ^E(^Hj dt = J tf r l|v(x(t),f)||ctt, (2) 

where (x\(t), X2(t), ...,x n (t)) denote the components of the trajectory x(t) in W 1 . Clearly, 
Mi depends on the initial point x* and the time interval [t* — r, t* + r] (and the vector 
field v). 

A Heuristic Argument. At this point it is insightful to give a heuristic argument as 
to why Mi is useful for revealing the geometric structures in the phase space of Q. Mi 
measures the arc-length of trajectories on a time interval (t* — r, t* + r). Trajectories 
with "close" initial conditions that remain "close" as they evolve on the time interval 
(t*—T, t* + r) are expected to have arc-lengths that are "close" . However, at the boundaries 
between regions comprising trajectories with qualitatively different behaviour in the course 
of their time evolution over the time interval (t* — r, t* + r), we expect the arc- lengths of 
trajectories starting on either side of the boundary to not be "close" on the time interval. 
(Of course, this will depend on r, and the dependence on r is an issue that we will shortly 
address.) Hence, these boundaries are denoted by an "abrupt change" in Ml, where "abrupt 
change" means that the derivative of Mi transverse to these boundaries is discontinuous 



on the boundaries. This will be more precisely discussed and defined in Section 2.1 Such 
boundaries correspond to the stable and unstable manifolds of hyperbolic trajectories, if 
we are in a hyperbolic region, of the phase space. Hence, in this way Mi detects stable 
and unstable manifolds of hyperbolic trajectories. 
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A General Construction for Lagrangian Descriptors. The arc length of a segment 
of a trajectory, defined by Q, is obtained by integrating the modulus of the velocity 
(||v||) along a trajectory for a time that parametrizes the curve. However, consistent with 
the heuristic argument given above, it is reasonable to expect this methodology to be 
successful for other positive scalar valued functions that embody an intrinsic physical or 
geometrical property of trajectories that are integrated along trajectories over the time 
interval (t* — r,t* + r). For example, integration along trajectories of other positive scalar 
valued quantities could have been considered, such as the modulus of acceleration (| |a| |), the 
modulus of the time derivative of acceleration (||da/dt||), or the modulus of combinations 
of v, a or da./dt, with the only restriction being that the integrals of these quantities along 
trajectories are well-defined. The heuristic argument would still apply and indicate that at 
the boundaries of regions comprising trajectories with qualitatively different time evolutions 
the accumulated value of the chosen positive quantity will "change abruptly" in the sense 
that the derivative transverse to the boundary of the corresponding Lagrangian descriptor 
will be discontinuous on the boundary, and examples of such regions would be those that 
are separated by the stable and unstable manifolds of hyperbolic trajectories. For this 
reason abrupt changes in a Lagrangian descriptor, which are manifested as discontinuities 
on the boundary of the derivative of the Lagrangian descriptor transverse to the boundary, 
are seen as singular curves in contour plots of the Lagrangian descriptor, and these are 
expected to be related to invariant manifolds. We will see specific examples that make 
these ideas precise in the remainder of this section. 

However, for now we note that based on this general reasoning a method for constructing 
families of Lagrangian descriptors for general time dependent flows would be as follows. 
Let | denote the bounded, positive intrinsic physical or geometrical property of the 

velocity field that is of interest. Consider a trajectory x(t) satisfying x(i*) = x* and defined 
on the time interval (t* — r,t* + r) (where r is "appropriately chosen" , as we will shortly 
discuss. Then |J r (x(t))| is a scalar valued function of t that depends parametrically on x* 
and t* . Therefore for 7 < 1 we can consider its L 7 norm: 

M(x*,r) Vir = (£ J^(x(t))PdfcV (3) 
Alternatively for 7 > 1 the L 7 norm is given by: 

M(x*,t* ) VtT = ^jf ^ FMtwA 7 • (4) 

For example, for ^ we considered the L l norm of the magnitude of the velocity evaluated 
on trajectories defined on the time interval (t* — r,t* + r). Other choices are possible. 
For example, we could consider the L 7 norm of the magnitude of the acceleration. For 
this case, for 7 = 1, we denote the resulting M function by M2(x*, t*) v>T . We could also 
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consider the L 7 norm of the magnitude of the velocity or the acceleration. In this case we 
denote the resulting M function by M^(x* , t*) VjT , and we will consider M3 for both 7 = 2 
and 7 = 5- Another possibility is to consider the L 1 norm of the magnitude of the time 
derivative of the acceleration, and in this case we denote the M function by M^x* , t*) v>T . 
We will explore more deeply the advantages and disadvantages of these different choices 
in the remainder of this section. However, we remark that the L 7 norms, 7 > 1, are not 
useful as Lagrangian descriptors in the sense that they do not reveal significant underlying 



structures in the flow. We will explore the reasons for this in Section 2.1.1 



Finally, the curvature is an intrinsic property of curves. It is a positive quantity that 



combines both velocity and acceleration (Kreyszig (1991)): 



_ V(w)(a-a) - (va) 2 

(vv)3/2 [b) 

The value of the curvature ranges from zero to infinity-the curvature of straight lines is 
zero and the curvature of fixed points is zero. We define a Lagrangian descriptor based on 
curvature as follows: 

|/i(x(t))l = (6) 

where singularities are avoided by introducing a > in the denominator. We have nu- 
merically verified that good choices for 'a' lie in the interval 5 > a > 1. This interval is 
determined considering that if 'a' is very large then Eq. Q tends to be independent of 
K,(x(t)). On the other hand if 'a' is very small the expression ([6]) will present very sharp 
profiles for zero curvatures near straight lines. In the examples to follow we will present 
results for a = 1. 

A key property of all of the Lagrangian descriptors is that they are quantities that 
"accumulate" along a trajectory, i.e. they are integrals of a positive quantity along a 
trajectory. The heuristic argument would not be valid for the integral of a quantity that 



changes sign along a trajectory. We will address this issue explicitly in Section 2.3. In Table 
[2] we summarize the different Lagrangian descriptors that we have introduced in this section 
and Appendix A describes in detail issues associated with the numerical computation of 
Lagrangian descriptors. 



Hyperbolic and Elliptic Regions. The phrases "hyperbolic regions" and "elliptic re- 
gions" occur throughout the literature. Here we want to consider their meaning more 
carefully. 

Traditionally, the words "elliptic" and "hyperbolic" have referred to the stability char- 
acteristics of individual trajectories. Trajectories are said to be hyperbolic if none of its 
Lyapunov exponents are zero (with the exception of the Lyapunov exponent tangent to the 
direction of motion of the trajectory, i.e. the Lyapunov exponent in the direction of the 
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Lagrangian Descriptor 


Integrand 


Norm 


Mi 


magnitude of velocity 


L 1 


M 2 


magnitude of acceleration 


L 1 


M 3 


magnitude of acceleration or velocity 


i 9 
La or L z 


Mi 


magnitude of the time derivative of the acceleration 


L 1 


M 5 


positive quantity related to curvature ( Eq. ( 6 ) ) 


L 1 



Table 1: The left hand column denotes the Lagrangian descriptors used, the middle column 
describes the positive quantities used as the integrand for the Lagrangian descriptors, the 
right hand column describes the norm that is used. 



velocity). Exponential dichotomies are also, equivalently, used to characterise hyperbolic- 



ity, see Coppel (1978); Dieci et al. (1997); Dieci and Vleck (2002). These characterisations 



of hyperbolicity are independent of the nature of the time dependence. The definition of 
"elliptic" for general time dependent trajectories has not received the same level of atten- 
tion, at least in the applied areas. Certainly, the notion of "elliptic stability type" for time 
periodic trajectories is well established. A time periodic trajectory is elliptic if all of its 
Floquet multipliers lie on the unit circle (there are some issues when a multiplier is 1 or 
— 1, but that level of technicality is not important for our discussion). 

However, there is a significant body of literature that does address the notion of "elliptic 
stability" for trajectories in vector fields having arbitrary (e.g. aperiodic) time dependence. 
It would be very interesting to explore the implications of these ideas in applications, 
in much the same way that Lyapunov exponents have been explored and exploited in 
applications. However, here we are considering what is meant by the phrases "hyperbolic 
region" and "elliptic region". 

With the notion of hyperbolic and elliptic stability of individual trajectories in hand, 
we can now address the question of what is meant by the phrases "hyperbolic region" and 
"elliptic region". Initially, one might guess that these are regions (e.g. open sets) having 
the property that all trajectories in such regions are either elliptic or hyperbolic. How- 
ever, there is an issue with this particular definition that requires careful consideration. 
Namely, traditionally notions of "elliptic" or "hyperbolic" stability type are infinite time 
characterisations of trajectories. Hence, in order to establish that all the trajectories in a 
region are of a particular stability type, one would need the trajectories to remain in that 
region for all time. In other words, one would need the region to be invariant (at least in 
forward time) . This requirement could (naturally) be met for linear or integrable Hamilto- 
nian systems, but not for general dynamical systems. There is an additional complication 
that occurs for "generic" dynamical systems (at least for dynamical systems defined by 
time periodic velocity fields). General theoretical results imply that it is unreasonable to 
assume that we can "cleanly" separate dynamical behaviour into "hyperbolic" and "ellip- 
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tic" . In particular, conservative generalisations of the "Newhouse phenomena" Newhouse 



(1977 


); 


Duarte 


(1999); 



stable and unstable manifolds of hyperbolic trajectories are "persistent" and will give rise 
to elliptic "islands of stability". If one "perturbs away" such tangencies, then new ones 
will be created elsewhere. We emphasis that, to date, such results have only been proven 
in the time-periodic setting. 

Despite ambiguities, it has nevertheless proved useful in the fluid mechanical situation 
to attempt to partition the flow into "elliptic" and "hyperbolic" regions. Early attempts in 
the fluids community to achieve such a decomposition of the flow led to the Okubo- Weiss 



criterion (Okubo (1970); Weiss (1991)), which is essentially Eulerian in nature. Efforts to 



generalize the idea to a Lagrangian setting appeared in Haller (2001) (but see also Due 



and Siegmund (2008); Branicki and Wiggins (2010)). This is not a problem that we will 
explicitly address in this paper. Nevertheless, we will use the phrases "hyperbolic region" 
and "elliptic region", but for most usages we will be considering linear and integrable 
systems, where the difficulties mentioned above are not an issue. 



2.1 Lagrangian Descriptors in Hyperbolic Regions 

Hyperbolic regions contain the "seeds" of change, uncertainty, and chaos in flows. In 
particular, they contain hyperbolic trajectories and their stable and unstable manifolds. 
In this section we will explore the performance of the Lagrangian descriptors defined in the 
previous section in familiar examples that exhibit hyperbolic behaviour. We will also relate 
their performance to more familiar diagnostics for "diagnosing" phase space structure, such 
as finite time Lyapunov exponents (FTLEs) and the ergodic decomposition. 

2.1.1 The Linear Saddle Point 

The first example that we analyze is the linear saddle. Initially, one might consider this 
example as "too trivial". However, we will show that it yields some interesting insights. 
Moreover, it lends itself to either exact, or accurate approximate, solutions for several 



interesting quantities in Section 2.3 



The velocity field that we consider is given by: 



x 

y 



Ax, 

-Ay, A > 0, 



(7) 



and the flow generated by this velocity field is given by: 



x(t,x ) 

y(t,yo) 



xoe 

yoe 



-xt 



A > 0, 



(8) 
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ll 1 
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Figure 1: Evolution of the contours as r increases, a) r = 0.5, b) r = 2, c) r = 10, d) 
isolated contours showing the evolution for r = 0.5, 2, 10. 

The origin, (x, y) = (0, 0) is a hyperbolic fixed point with stable and unstable manifolds 
given by: 



VT(0,0) = {(x,y) eM 2 |x = 0, y^0}, (9) 
W u (0,Q) = {(x,y) £ R 2 \y = 0, x / 0}. (10) 

In Figure [T] we show the evolution of the contours of M\ for r = 0.5,2, 10. The first 
thing to observe is that the patterns of the contours displayed depend on r. For r small 
the structure of M\ is smooth and for increasing r the contour patterns converge towards 
a structure that displays the manifolds by means of discontinuity in the derivatives. An 



analytical argument concerning the convergence time is given in Section 2.1.2 

In Figure [2] we show several different Lagrangian descriptors, where each is computed 
for r = 10. Figure ^a) shows contours of the Lagrangian descriptor M\ for A = 1. We see 



that the stable and unstable manifolds of the hyperbolic fixed point at the origin are easily 
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Figure 2: Contour plots for a selection of Lagrangian descriptors for ([7]). a) Mi; b) M3, 
using the La norm; c) M3 using the L 2 norm; d) M5. 



identified. For A = 1 it is easy to see from the simple form of the velocity field ([7]) and its 
flow given in Q that M\ = M2 = M4, and therefore the latter two Lagrangian descriptors 
are not shown. 

Figure [2^) shows the contours of M3 for 7 = 0.5-the Ls norm of the modulus of the 
acceleration-for which the stable and unstable manifolds of the hyperbolic fixed point at 
the origin are readily apparent. However, figure [2)3) shows the contours of M3 for 7 = 2- 
the I? norm of the modulus of the acceleration, and in this case the invariant manifolds 



we 



structure is completely lost. We explain this behaviour in Section 2.1.2. In Figure [2 
show contours of the Lagrangian descriptor M5 (for a = 1). In this case the manifolds are 
clearly shown. This is not surprising since the stable and unstable manifolds of the origin 
are the only trajectories having zero curvature. In fact, in this case we would expect that 
the contours of M5 would reveal the manifolds for relatively small r. In particular, we have 
verified numerically that essentially the same contours are obtained for r = 2. 

Finally, we note that the Lyapunov exponents can be computed analytically for every 
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trajectory in this example. Since the system is linear they are the same for every trajectory- 
±A (except for the exception case of trajectories starting on either the stable or unstable 
manifold of the origin, in which case one of the Lyapunov exponents will be zero). Hence 
the contours of constant Lyapunov exponents reveal no structure in this case (as noted in 



Branicki and Wiggins (2010)). This will be explored in more detail in Section [273 



2.1.2 Analytical Justification 

In this section, for the specific example ([7]), we give analytical justification underlying the 
heuristic argument as to why contours of the Lagrangian descriptors having the property 
that the derivative of the Lagrangian descriptor transverse to the contours correspond to 
material curves. This also provides insight as to why the IS 1 norms, 7 > 1, do not reveal 
this singular property in the contours for Lagrangian descriptors that use such norms. 
For ([7]), using ([2]), the M\ function is given by: 

Mi(xo,yo;*o,r) = y/x(t,x ) 2 + y(t,y )*dt = f 0+T \/x 2 X 2 e 2Xt + y 2 X 2 e~ 2Xt dt (11) 

This function is clearly nonnegative since it is a measure of arc length of a trajectory. It is 
clear by inspecting ^ that M\(0, 0; to, r) = 0. This makes sense since the arc length of a 
point is zero. We want to compute the integral M\. For simplicity we assume without loss 
of generality that to = (this is always possible for an autonomous system after a time 
shift t' = t — to and by redefining the initial conditions xq and yo). 

As already noted the patterns of the contours displayed in Figure [I] depend on r. For 
r small the contours are smooth but for increasing r they develop singular features along 
the unstable and stable manifolds. The nature of what we mean by the phrase "singular 
features" will be explained in the following where we provide analytical justification that 
underlies this behaviour. 



Behavior for Small r. We now want to argue that for r small, (11) does not display 
"sharp features" that are related to the stable and unstable manifolds of the hyperbolic 
fixed point of ([7]). Towards this end, we consider an expansion of the expression (11) for 



fixed yo that is valid for xq small {i.e. near the stable manifold) and for small r. The 
expansion has the following form: 



Mi(aro,l/o;io,r) = f ^ (x X) 2 e 2 ^ + (y X) 2 e- 2Xt dt 

= f^X\y \ + ^^x^dt + O(xt). (12) 
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Figure 3: In blue, the numerical evaluation of M\ along xq for yo = 0.5, r = 2, A = 1. In 



red, the evaluation of the leading order terms (13). 



It is important to notice that in the integrand xq is always multiplied by e 2Xt , thus if we 
want the "error terms" in the series expansion around xo = to be small, then e 2Xt should 
not be large, which means that r should not be large, as we are assuming. The leading 



order terms in (12) can be explicitly computed, and are found to be: 



M lt i o (x ,y ;t ,T) 



,-At 



AM + 



3At 



Mvo\ 



-At 



\yo\ + 



ii | dt, 

~ 3At ) \vo\ 



(13) 



For yo fixed, (13) is a smooth function of xq. Figure [3] shows a comparison of the accuracy 
of (12) with (13) for r = 2 along a horizontal line perpendicular to the stable manifold. 



More precisely, we fix i/q = 0.5 and plot (12) (computed numerically, and shown in blue) 



and (13) (shown in red). It is clear from the figure that leading order terms (13) accurately 



approximate (|12|) in the interval xq £ [—0.01,0.01]. 



this end, we re-write (11) as follows: 



Behavior for Large r. Now we consider the behaviour of ( 11 ) when r is large. Towards 

(14) 



Mi(x ,yo;t ,T) 



X 2 + Y 2 dt + 



X 2 + Ygdt 



where we have defined Xq = XgA 2 e 2A * and Y 2 = yl\ 2 e 2X1 
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In the first integral term in (14), since t is positive if r is large Xq may assume large 
values for t values approaching the upper limit of the integral. By similar reasoning, Yq 
will become small for t values approaching the upper limit of integration. For the second 



integral term in (14), since the integration range is negative Yq assumes large values for t 



approaching the lower limit of the integral. By similar reasoning, Xq will become small for 
t values approaching the lower limit of integration. These facts suggest isolating the large 
parts of both integral terms by choosing constants t p , r n > (as we will see, precise values 



for these constants are not important for the argument) and re-writing (14) as follows: 



Mi(x ,yo;io,r) 



+ 



X 2 + Y 2 dt + 



X 2 + Y 2 dt + 



X 2 + Y 2 dt 



X 2 + Y 2 dt 



(15) 



We expand the first term in (|15|) in and the second term in (15) in and rearrange 



the terms to obtain the following: 



Afi(x ,W);to,T) =/ \X \dt+ \Y \dt + O(l/X ) + O(l/Y ) + B, (16) 

J Tp J —T 

where B is defined as: 



B 



X 2 + Y 2 dt + 



X 2 + Y 2 dt. 



(17) 



Using the fact that for r large |xo|e Ar >> |xo|e ATp and |yo|e AT >> |yo|e Ar " it is straightfor 



„At 



,At„ 



ward to determine that the leading order terms of ( 16 ) are given by: 



Mi } i (x ,yo;t , r) = |x |e T + |y |e 



(18) 



In the above expressions it is important to note that the estimates are valid even if xq 
or yo are small as far as the products XQe xt , yoe~ xt are large. This occurs very close to the 
origin, xq ~ or yo ~ 0, if r is large, and this corresponds to a large integration time in 
the definition of the Lagrangian descriptor. 

In this case discontinuities in the derivatives (or even on the functions themselves in 
other possible examples) observed along the stable and unstable manifolds rise because the 
asymptotic expansions at both sides of the manifolds, do not match well across these lines. 



In practice for any finite r the series expansion (12) is still valid, thus strictly speaking the 
function M is regular, but it is accurate only in an extremely narrow gap around xo = 0, 
exponentially decreasing for increasing r. 

Figure [4] shows a comparison of the accuracy of (14) with (18) for r = 10 along a 
horizontal line perpendicular to the stable manifold. More precisely, we fix yo = 0.5, 
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Figure 4: In blue the numerical evaluation of (14) along xq for yo = 0.5, r = 10, A = 1 



In red the evaluation of expression (18). 



t = 10, A = 1 and plot (14) (computed numerically, and shown in blue) and (18) (shown in 



red). It is clear from the figure that leading order terms ( 18 ) accurately approximate ( 14 ) in 
the interval xq £ [—0.1, 0.1] and that the derivative of both functions appears discontinuous 
on the stable manifold. 



The L 7 norm, 7 > 1. In Section[2]we stated that the contours of the L 7 norm, 7 > 1, of a 
bounded, positive intrinsic physical or geometrical property of the velocity field integrated 
along trajectories do not yield interesting structures in the flow. Here we give an argument 
why this is the case in the context of the example Q. 

We first extend the previous results for the L 7 norm, with 7 < 1. The essential feature 
to consider is the exponent on the integrand, which we will take as the magnitude of 
velocity for this example: 



M(x ,y ;t ,T) = 




Following the same argument as above, in this case the asymptotic expansion is: 

M(x , y ; t , r) ~ J* (x 25 + 5Y X 2 { - 1+5) ) dt + J ^ (y 2S + 5X Y 2{ - 1+S) ) dt 

+ o(i/x 4 - 25 ) + o(i/y 4 - 25 ) (20) 

Following an argument similar to that given above, for < 5 < 1 the leading order terms 
in the expansion are: 
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M (X , 2/o; t0, T) ~ M e 7Ar + JyoT g7 Ar (21) 

7 7 

where 7 = 25. If 7 < 1 the derivative across the stable manifold is singular, which explains 
the enhancement of the features displayed in Figure [2)3) for this choice of exponent. The 
case 7 = 1 corresponds to Q. On the other hand, if 7 > 1 the expression above should be 
raised to the power I/7. In this case the derivative across the stable manifold is continuous. 
We see this in Figure [2p) (7 = 2) for which the contours of M are smooth. 



2.2 Lagrangian Descriptors in Elliptic Regions 

Now we will consider the behaviour of Lagrangian descriptors in "elliptic regions" . In the 



same spirit as in the discussion in Section 2.1.1, we will begin by discussing the familiar 
linear elliptic fixed point. 

2.2.1 The Linear Elliptic Point 

We consider the velocity field: 

x = y, 

V = -x, (22) 
and the flow generated by this velocity field is given by: 

x(t,xo) = xq cost + yo sin t, 

y(t,yo) = -xo sin t + y cost. (23) 

The origin, (x, y) = (0,0) is an elliptic fixed point. In this case the M\ function has the 
form: 

rto+r 

Mi(x ,y ;to,T) = / ^x{t,xo) 2 + y(t,y ) 2 dt 

J tQ-T 

pto+r 

= / \J (— xo sin t + yo cos t) 2 + (— xq cos t — yo sin t) 2 dt, (24) 

J to—T 

which is easily computed to give: 

M 1 (x ,ya;t ,T)= f ^x 2 + y 2 dt = 2T^x 2 + y 2 . (25) 

which is a smooth function for all r except at the origin. The contours of Mi, for any fixed 
r, are smooth circles surrounding the origin. Also, note that for this example it is easily 
verified that M t = M 2 . 
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Figure 5: Evaluation of Lagrangian Descriptors for the system (27). a) M± for r 
h)M\ for r = 75 ;c) M3, using |a| in the integrand with the norm. 



25; 



2.2.2 An Integrable Example: Shear and Resonance 

In this section we consider an example that, in some sense, is the "simplest" nonlinear 
problem that has entirely elliptic behaviour-a one degree-of-freedom integrable Hamil- 
tonian system in action-angle variables. Nevertheless, this system provides some useful 
insights that raises questions for more complex systems. 
The Hamiltonian is given by: 

H(I) = T - -, (26) 
and the associated Hamiltonian vector field is given by: 



dH 

~dT 
dH 

~d6 



I 6 -I, 



0. 



(1,0) elxS 1 , 



(27) 



It is evident from (27) that I is constant, and each value of / corresponds to an invariant 



circle, or "invariant one torus" (since 6 is periodic). The frequency of each torus is given 
by ^j- = I 3 — I and the shear, or "twist" associated with each torus is ^f- = 3I 2 — 1. The 
typical behaviour for each / = constant invariant circle is that it has nonzero frequency 
and nonzero shear. However, there are five exceptional tori that we note. 

Resonance. The three invariant circles 1 = 0,1 = ±1 correspond to circles of fixed points 
since the frequency is identically zero on each of these circles. 

Zero shear. On the invariant circles / = ±i the frequencies are -F^tj but the shear, or 
"twist" is identically zero on each of these invariant circles. In the dynamical systems 
literature these are referred to as "twistless" invariant circles. The choice of "twist 



less" tori will allow us to compare with recent work in Beron-Vera et al. (2010) when 
we consider FTLEs in Section [2.31 
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Figure [5p) illustrates how Mi detects the three circles of fixed points (dark blue) as 
well as the high shear "twistless tori" (the reddish colours), for r = 25. Figure]^) shows 
the contours of M for r = 75, which contain essentially the same information as for r = 25. 

The function M3, based on |a| and the L2 , gives a different perspective on the structure 
of the flow. Since in this system all trajectories move with constant velocity, their accel- 
eration is zero, from which it follows that M3 is identically zero. In this case Figure |5p) 
confirms a flat structure which is obtained for all r > 0. That is, there is no convergence 



threshold in r. We will explore this example again when we discuss FTLEs in Section 2.3 



2.3 Remarks on the Application of Finite Time Lyapunov Exponents 
(FTLEs) and Time Averages for These Examples 

In this section we want to compare two other techniques that are used to visualise phase 
space structures-finite time Lyapunov exponents (FTLEs) and the ergodic decomposition- 
with the method of Lagrangian descriptors. 



2.3.1 Finite Time Lyapunov Exponents (FTLEs) 



FTLEs have proven to be a useful tool for obtaining a "global picture" of phase space struc- 
ture. The fundamental theorem concerning Lyapunov exponents (i.e. for infinite time) is 



the Oseledets multiplicative ergodic theorem (Oseledets (1968)). This theorem gives suffi- 
cient conditions for finite time Lyapunov exponents to converge to their infinite time values. 
However, it does not provide information how accurately FTLEs, for a specified finite time, 
approximate their asymptotic values, nor does it relate Lyapunov exponents (either finite 
or infinite time) to invariant manifolds (e.g. "material" curves for two-dimensional, time 
dependent velocity fields). Discussions of Oseledets theorem in the context of applications 



can be found in Legras and Vautard (1996); Lapeyre (2002). 



FTLEs have been used successfully as "proxies" for invariant manifolds ( Haller ( 2000 ) ; 
Shadden et al. (2005)), in the following sense. An initial time is chosen and the spatial 



domain is decomposed into a discrete grid. These grid points are then integrated forward 
(resp., backward) in time, for a fixed time (each grid point is integrated for the same 
time). The maximal Lyapunov exponent for each trajectory is computed for this time (the 
"FTLE") and a value is given to the initial point corresponding to this maximal FTLE. 
This is done for each grid point and the result is a forward (resp., backward) FTLE field 
at the chosen initial time. Then the level curves of the FTLE are displayed-different 
colours for different values. Locally maximal level curves of the forward (reps, backward) 
FTLE field are typically identified with the unstable (resp., stable) manifolds of hyperbolic 
trajectories. In the "topographical map" of the level curves of the FTLE field the locally 
maximal level curves have the property of a "ridge", i.e., moving in a direction that is 
not tangent to the locally maximal level curve results in moving in regions of smaller 



FTLEs. The quantification of this notion of a "ridge curve" is carried out in Shadden et al. 
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(2005), and we mention this again below. For now we note that the identification of locally 



maximal level curves of the FTLE field has been noted to provide good approximations 
to material curves in certain examples. Nevertheless, there are numerous issues , the full 
details of which are sometimes not reported in papers, that use FTLEs for specific physical 
applications. In particular, we point out three common issues. 

The choice of finite time. This is a problematic issue for which there is no theoretical 
guidance. The Oseledets theorem gives sufficient conditions for the FTLE to converge. 
However, for any finite time, it is not clear how "close" the FTLEs computed at that 
time are to their asymptotic values. The limited theoretical results that are available 
predict a "relatively slow" convergence time (e.g. logarithmic in r, in our notation, 
see, e.g., Goldhirsch et al. (1987)). Moreover, it is observed in the course of time 



evolution trajectories that are "finite time hyperbolic" or "finite time elliptic" can 
change their "finite time stability type" (see, e.g. Branicki and Wiggins (2010); 



Mancho et al\ ( |2006bD ). There are numerous issues related to this topic that demand 



further theoretical investigations. 

Discretization issues. Once the FTLE fields are computed the goal is then to look at 
the level curves of the FTLE field. However, for any realistic system, FTLEs must be 
computed numerically, requiring discretization in space (and possibly time). Setting 
the issue of the choice of finite time aside, the computation of FTLE fields may be 
somewhat noisy, and this means that the "raw output" of FTLE computations does 
not yield smooth level curves. In general, some form of filtering and smoothing is 
required in order to obtain a FTLE field that yields smooth level curves. This can be 
a subjective process. Examples of the effects this can yield, in the context of known 



"benchmark" examples can be found in Branicki and Wiggins (2010). 



The connection to invariant manifolds ("material curves"). Setting aside the cru- 
cial issue of the smoothness of the level curves of the FTLE field, it is often stated 
in the literature that the level curves that locally maximal level curves of the FTLE 
field are "material curves" . In general, this is not true, although they may be "close" 
to material curves and provide a good approximation for methods that compute ma- 



terial curves via particle advection. This issue was carefully considered in Shadden 



et al. (2005), who introduced the notion of a ridge curve in the FTLE field as an 



approximation to a material curve. They showed that ridge curves tend to be a good 
approximation to a material curve by deriving an expression for the flux across the 
ridge curve, which may be small, but generally nonzero. There is a great deal of 
ambiguity in the literature concerning the notion of a "ridge curve" or "ridges of the 
FTLE field" . In many papers the notion of a "ridge" is taken to by synonymous with 
locally maximal level curves of the FTLE field. However, this is not the case, and 



ridge curves are precisely defined in Shadden et al. (2005). 
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Of course, a goal of this article is not to present a critical discussion of FTLEs and the 
range of their applicability (this is dealt with, to varying degrees, in Mancho et al. (2006c); 



Branicki and Wiggins (2010)). Rather, we will compute FTLEs, for different times, for the 



simple examples that we have previously considered and compare them with the results 
obtained from Lagrangian descriptors. While these examples are certainly far simpler than 
a typical "realistic" geophysical fluid dynamics application, they have the advantage that 
the FTLEs, hyperbolic trajectory (if it exists) and their stable and unstable manifolds are 
known analytically and, in that sense, they may serve as benchmarks for testing various 
methods. 

For linear, time independent velocity fields the (infinite time) Lyapunov exponents are 
simple to compute. They are the real parts of the eigenvalues of the matrix associated 
with the linear velocity field. Therefore, for the linear saddle point, Q, the Lyapunov 
exponents are ±A (with the maximal Lyapunov exponent being A > 0) and for the linear 
elliptic point, (22), the Lyapunov exponents are both zero. The FTLE fields for each of 
these examples is shown is Figure [6] for r = 10, 25, 75. For each r a "flat field" is shown, 
which is to be expected for linear velocity fields since the Lyapunov exponents are the same 
for all trajectories. Hence, the FTLEs reveal no structure in the flows, in contrast to the 
Lagrangian descriptors. 

In Figure [7] we compute the FTLEs for (27) for r = 25 and r = 75. The velocity 
field is integrable and expressed in action-angle coordinates. Therefore we know that the 
Lyapunov exponents are identically zero. However, since the velocity field is nonlinear, the 
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Figure 7: FTLE for (27). a) r = 25; b) r = 75. 



approach to the limit is "nonuniform" in the sense that different trajectories show approach 
the limit at different rates. Consequently, Figure [7] shows "structure" in the FTLE fields 
(in contrast to our example linear velocity fields), but this structure will disappear in the 
infinite time limit. It is clear from the colour bar in Figure [7] that the FTLE are decreasing 
with r. Note that the smallest values for the FTLE appear near the invariant circle 
I -- 



n 




which correspond to the "twistless" invariant tori. It was shown in 



Beron-Vera 



et al. (2010) that twist less tori in integrable systems expressed in action-angle variables 



]0.06 
0.04 
0.02 



will have zero FTLE, for any time. 



2.3.2 Time Averages 



In this section we will describe another approach for discovering phase space structure 
in dynamical systems that utilises time averages of certain functions along trajectories. 
Relating time averages to phase space structure is reminiscent of notions of ergodicity, 
and a general framework that has been used in the context of fluid flows (among other 



applications) is the ergodic decomposition. This approach was developed in Malhotra et al. 



(1998); Poje et al. (1999); Mezic and Wiggins (1999), and is based on fundamental work of 



Rokhlin (1966). Superficially, this approach appears to have some similarity to the method 
of Lagrangian descriptors. We will explain that this is not the case, and consider the 
application for the linear saddle Q, where this approach fails in several aspects. 

We begin by giving a very brief description of the method. The setting is that of 
smooth ergodic theory. We have a compact domain, A, and a differentiable dynamical 
system defined on A possessing an invariant measure. The dynamical system can be either 
discrete time, i.e. a map, or continuous time, i.e. a flow, and by "flow" we mean the precise 
mathematical definition as a one-parameter group of transformations of A. We denote by 
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L 1 (^4) the space of scalar valued functions on A having the property that the integral of its 
absolute value over A exists. It is essential that the Birkhoff ergodic theorem is applicable. 
In particular, we require that the time average of functions / G L l (A) along trajectories 
of the dynamical system exist for almost all (with respect to the invariant measure) initial 
conditions in A, and we point out that by "time average" we mean infinite time average 
(as in the statement of the Birkhoff ergodic theorem). 

Suppose we choose / G L l (A), and consider the function on A defined by the time 
average of / along all possible trajectories in A. It can be shown that the level sets of this 
function are invariant sets with respect to the dynamics. Now suppose we have an inner 
product on A and we are able to choose a set of mutually orthogonal functions, {/«}, iEN, 
having the property that the set of all possible linear combinations of the {/j} are dense in 
L 1 (^4). For each /j, we consider the function of A given by the time average of trajectories 
through initial points in A, which we denote by /*. Finally we consider the "joint level 
sets" that are constructed by considering the intersections of all level sets of all of the /* . 
This collection of "joint level sets" defines a partition of A, and the dynamical of each 
element of the partition is ergodic (the "ergodic partition" ) . 

An essential tool used in this procedure is the Birkhoff ergodic theorem which requires 
compactness of the phase space as an essential ingredient. Unfortunately, this theorem 
has not been proven for time dependent systems with general, aperiodic, time dependence. 
The theorem holds for maps and for time periodic and quasi periodic velocity fields. In 
the latter case the dimension of the system can be increased by considering time as an ad- 
ditional dependent variable in the definition of the velocity field. Since the velocity field is 
quasi periodic (of which periodicity is a special case) time is treated by including the inde- 
pendent (but finite) angles associated with each frequency component of the quasi periodic 
time dependence. In this way compactness is preserved (each included angular variable 
is compact) and the velocity field, in this higher dimensional setting, is autonomous and 
generates a flow on a compact space. For general, aperiodic time dependence transform- 
ing the problem to an autonomous setting in which compactness is preserved is generally 
not possible (see, e.g. Balibrea et al. (2010) for a recent review of the relevant issues for 
aperiodic time dependence). 

Even though the Birkhoff ergodic theorem has not been proven for general aperiodic 
time dependence, the general approach suggested by the ergodic decomposition has proven 
fruitful, especially when the /, under consideration have a physical meaning (and this is 
not unlike the approach for constructing Lagrangian descriptors). In particular, the finite 
time average of the horizontal component of the velocity field has been shown to provide 



insight in the flow in numerous different situations by Mezic et al. (2010); Poje et al. 
fll999D ; |Malhotra et al.\ fll998| ). For the general velocity field, (|lj), let x(i, t*,r) denote the 
trajectory satisfying x(x(i*, £*, r) = xo and let v x denote the x component of the velocity 
field. Then the finite time average of v x along this trajectory is denoted by: 
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v*(xo,i*,T) = 1 ^ +T v x (x(t,t*,T),t)dt (28) 

T Jt* 

Since the integral is only over a finite length (in time) of the trajectory, conditions for its 
existence are relatively mild. However, the "meaning" of such results for such finite time 
averages is not a priori clear. The lack of a rigorous theoretical framework is similar to 
the situation for FTLEs. The typical approach is to compute the relevant quantities (e.g., 



(28)), "see" if the computations yield something "interesting", and then to investigate 
further. 

The linear velocity field ([7]) is defined on a non compact domain, thus the Birkhoff 
ergodic theorem does not apply to it, in the sense that it does not allow us to conclude 
that time averages along trajectories exist. Nevertheless, we can still compute the finite 
time averages along trajectories of physically interesting quantities in the spirit of the 
patchiness work described above. For this example the finite time average of the horizontal 
component of the velocity field can be explicitly computed. Taking t* = 0, A = 1, and 
denoting the finite time average of v x by v*, we have: 

v * = ^exp(r)-l). (29) 

Clearly, the limit does not exist as r — > oo. In figure [8] we plot the contours of (29) for r = 
10. It is also clear that the contour plot does not reveal the presence of the stable manifold 
nor invariant sets, despite the system, which is hamiltonian, contains invariant sets defined 
by H (x, y) = Xxy = c . The stable manifold of the hyperbolic trajectory at the origin is 



given by xq = 0. We see from (29) that the finite time average changes sign at xq = (and 



this is true for any r). However, the derivative of (29) with respect to xq is smooth across 
xq = 0. Hence, the finite time average of the horizontal velocity component has no singular 
features that would indicate the presence of material curves. A reason for this behaviour 
is that general finite time averages do not require that the integrand is non-negative along 
trajectories, which is an essential requirement for Lagrangian descriptors, and motivation 
for this was given in the heuristic argument. Typically integrals of oscillating quantities 
along trajectories lead to a distorted output which makes interpretation of the results, in 



terms for phase space structure, difficult. In certain cases it was noted in Poje et al. (1999) 
that for increasing averaging time oscillations lead could lead to average velocities that 
approached zero everywhere, and as a consequence, in this limit any spatial structure is 
completely lost. 

2.4 Application to the Forced Duffing Oscillator 

In this section we will discuss further the ability of Lagrangian descriptors to highlight 
manifolds. Also we will compare and contrast the behavior of Lagrangian descriptors with 
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Figure 8: Integral of the horizontal component of the velocity along trajectories of the 
linear saddle point ([7]) for r = 10. 

FTLEs and time averages of velocity components along trajectories using a familiar, but 
much more complicated nonlinear system-the forced Duffing equation: 



x 

y 



y, 

x 



x 3 + ef(t), 



(30) 



We will consider three cases. 



Integrable case. In this case we take e = in (30). The phase space structure of the 
resulting system is well-known. The system has a hyperbolic fixed point at the origin 
connected by a symmetric pair of homoclinic orbits. The region inside and outside 
the homoclinic orbits consist of continuous families of periodic orbits, and correspond 
to elliptic regions. 



Periodically time dependent case. For this case we set e = 0.1 and f(t) = sin(i) (30). 
The system has a hyperbolic periodic trajectory near the origin, referred to as a 



Mancho et al. 


(2003 


); 


Madrid and 



Mancho (2009); Ide et al. (2002)). 



Aperiodically time dependent case. The aperiodically time dependent forcing, f(t), 
has the form shown in figure [9j which his obtained from a trajectory of the Duffing 
equation in a chaotic regime. We have chosen e = 0.15, which results in an amplitude 
of the time dependence that is, roughly, and order of magnitude smaller that the 
time periodic case considered above. In this case the system also has a hyperbolic 
trajectory near the origin. 
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Figure 9: The time series of the function f(t) corresponding to aperiodic time dependence 
for the Duffing Equation. 



Singular Features of the Contour Plots of Lagrangian Descriptors. In our study 
of the linear saddle in Section 2.1.1 we argued that Lagrangian descriptors had two impor- 
tant properties. 1) They depend on the time of integration r. 2) After a certain convergence 
time they display "singular contours" that correspond to the stable and unstable manifolds 
of hyperbolic trajectories. By "singular contours", we mean contours having the property 
that the derivative of the Lagrangian descriptor transverse to the contour is discontinuous. 
Here we will show that these two properties also hold for the cases of the forced Duffing 
equation that we have considered. In figure 10 we show the stable and unstable manifolds 
of the hyperbolic trajectory and contours of the function M\ for r = 2 and r = 10. For 
both cases we see that the contours of M\ converge to the stable and unstable manifolds. 

Figure 11 further illustrates how singular features of the contours of M\ correspond to 
stable and unstable manifolds of hyperbolic trajectories. The figure shows values of M\ 
along the vertical line shown in Figure 10 :) . We see that the value of M\ "change abruptly" 



at the locations of the manifolds. Discontinuities in the derivative of M\ quantify the notion 
of "abrupt change" of M\. The latter are illustrated in the dashed line in the figure. 

In figure 12 a) we show segments of the stable and unstable manifolds of the hyperbolic 
trajectory near the origin. In figure 12 b) we show the same manifolds overlaid with the 
contours of M\ computed for r = 10. The manifolds align well with the "singular contours" 
of Mi. 



Lagrangian Descriptors vs FTLE. In figure [13] we show the stable and unstable 
manifolds of the hyperbolic trajectory near the origin (left most column), the forward time 
FTLEs (middle column) and the contours of the Lagrangian descriptor M\ for the three 
cases, and all for r = 10. In all cases the "singular" contours of M\ approximate the 
manifolds much better than the contours of the forward FTLE, since although the forward 
FTLE only captures the stable manifold of the hyperbolic trajectory near the origin, it 
displays artifacts in the FTLE field that have no Lagrangian interpretation, and that are 
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Figure 10: For the integrable case: a) the stable and unstable manifolds of the hyperbolic 
fixed point; b) contours of M\ for r = 2; c) contours of M\ for r = 10. For the periodically 
forced Duffing equation: d) Segments of the stable and unstable manifolds of the hyperbolic 
trajectory near the origin ; c) contours of M\ for r = 2; c) contours of M\ for r = 10. 
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Figure 11: The black solid line represents the function M\ vs y at a fixed x. Mi along the 
vertical black line shown in Figure [Hfe ) . Abrupt changes in M\ , occurring at the location 
of the manifolds, correspond to discontinuities of the derivative of M\. The dashed line 
represents 0.1 times the derivative of M\ with respect to y. 
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Figure 12: The aperiodically time dependent case, a) The stable and unstable manifolds 
of the hyperbolic trajectory near the origin, b) An overlay of the stable and unstable 
manifolds of the hyperbolic trajectory near the origin with the contours of the M\. The 
manifolds and M\ are computed for r = 10. 



not present on the contour plots of M\. 



Lagrangian Descriptors in Elliptic Regions. Next we consider the behaviour of sev- 
eral Lagrangian descriptors in "elliptic regions" (cf. the discussion of elliptic and hyperbolic 
regions in Section [2]) for both the periodic and aperiodically time dependent cases. The 
elliptic region that we consider consists of the region enclosed by the left hand homoclinic 
orbit in the e = case (but we consider the behaviour of the Lagrangian descriptors in this 
region for the periodic and aperiodically time dependent cases). 

We first consider the Lagrangian descriptors M5 and M2 in the periodically time de- 
pendent case. Figure 14 a) and b) shows the contours of these Lagrangian descriptors 
computed for r = 70. The contours of each Lagrangian descriptor show a smooth struc- 
ture near the centre of the elliptic region. This is not surprising since for a "weak" time- 
periodic perturbation of an integrable system we expect "most" of the closed trajectories 
to be preserved as (two- frequency) Kolmogorov-Arnold-Moser tori. This would imply that 
trajectories starting "close" to each other would have a similar past and future on the time 
interval (t — r, t +t). Outside this smooth region the contours of M5 and Mi reveal a more 
complex structure. This "complex structure" appears to be associated with strong mixing. 
This suggests that r may be a useful parameter for quantifying the notion of "finite time 
mixing" , but making this notion precise requires further research. Practically, the notion of 
"finite time mixing" is very important in applications, and lies outside the scope of notions 
of mixing in traditional ergodic theory. 

Next we consider the aperiodically time dependent case in the same region. In this case 



27 




Figure 13: For the integrable Duffing equation: a) the stable and unstable manifolds of the 
hyperbolic fixed point; b) forward FTLE for r = 10; c) contours of Mi for r = 10. For the 
periodically forced Duffing equation: d) segments of the stable and unstable manifolds of 
the hyperbolic trajectory near the origin computed for r = 10 (and displayed at t = 0); e) 
forward FTLE for r = 10; f ) contours of M\ for r = 10. For the aperiodically forced Duffing 
equation: g) segments of the stable and unstable manifolds of the hyperbolic trajectory 
near the origin computed for r = 10 (and displayed at t = 0); h) forward FTLE for r = 10; 
i) contours of M\ for r = 10. 
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the isolated region of "trapped" trajectories shown in Figure 14 a) and b) is "broken" as it 



is "invaded" by segments of the stable and unstable manifolds of the hyperbolic trajectory 



near the origin, as shown in Figure 14 c) and d) for the Lagrangian descriptors M3 and 
M4 (more details are evident from the contours of M3). Although not shown, we note that 
the Lagrangian descriptors Mi, M2 and M5 give results similar to those of M4 and M3 in 
this case. 



Time Averages of Velocity Components. We now return to the issues raised in 



Section 2.3.2 In particular, we compute the finite time average of the horizontal component 



of the Duffing equation, following (28), for the three different time dependencies that we 
are considering. 

The results are shown in figure [l5j where the time averages are computed in each case 
for t = 10 and r = 70. For each time dependence we see similar behavior. For the smaller 
time we see "a hint" of structure that bears some resemblance to "fattened" versions of 
short segments of manifolds of the unstable manifold of the hyperbolic trajectory. For the 
longer time the contours of the finite time average of the horizontal component of velocity 
develop a more complex spatial structure that obscures the underlying unstable manifold 
structure of the hyperbolic trajectory. 

3 Applications to time dependent 3D flows 

In this section we show that Lagrangian descriptors can also provide accurate information 
on the stable and unstable manifolds of hyperbolic trajectories in three dimensional (3D) 
time dependent flows. Computation of the stable and unstable manifolds of hyperbolic 



trajectories in aperiodically time dependent flows was discussed in Branicki and Wiggins 



(2009), where an algorithm for their calculation was developed and several "benchmark" 
examples were considered. The particular example that we will consider is the perturbed 
Hill's spherical vortex, which we will take as a benchmark for the performance of our 
methods in 3D. We give a brief description of the velocity field. More details on the 



background of Hill's spherical vortex can be found in Branicki and Wiggins (2009). 
The velocity field, v, has the general form: 



v = H(x,y,z) + S(x,y,z,t) 



where H(x, y, z) is given by: 



H x = (u r sin + u® cos 0) cos (31) 
H y = (u r sin © + ue cos 0) sin <£, (32) 
H z = (u r cos — «e srn @)j (33) 
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a) 




Figure 14: Contour plots of the of Lagrangian descriptors in the region corresponding 
to that enclosed by the left hand homoclinic orbits for e = 0. Panels a) M5; b) Mi 
correspond to the periodically time dependent case. Panels c) M3; d) M4 correspond 
to the aperiodically time dependent case. The Lagrangian descriptors are computed for 
t = 70. 
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Figure 15: Contours of the finite time average of the horizontal component of the Duffing 
Equation along trajectories; The integrablegcase (e = 0) for a) r = 10; b) r = 70; the 
periodic case for c) r = 10 d) r = 70; the aperiodic case for e) r = 10; f) r = 70. 



Here: 

r = \/ x 2 + y 2 + z 2 , = arccos(z/r), <3? = arccos(a;/ \J x 2 + y 2 ) 



U(l - a 3 /r 3 ) cos 6, if r > a, _ J -C/(l + a 3 /(2r 3 )) sin G, if r > a, 
Ur ~ \ -f[/(l-r 2 /a 2 )cosG, if r < a, U@ ~ { p{l - 2r 2 /a 2 ) sin 6, if r < a, 



(34) 

The matrix S 1 is given by: 
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where 5 is the orthogonal matrix: 



cosip cos ip — sin tjj cos 8 sin cos ip sin + sin ip cos cos sin tp sin 
- sin tp cos — cos tp cos sin (p — sin -0 sin <p + cos ^ cos 9 cos 93 cos ^ sin 



sin c/ sin tp — sm 6/ cos ip 

and j4(t) is the time dependent function: 

A(t) = (0.3 + 0.27sin(3.3t)) exp(-(t - 6) 2 /2.5 2 ). 



Additionally: 



and 



= 0.5 + 0.05sin(2t), ip = 5t, ip = 



a (t) = -0.5, 0(t) = -0.5 and 7(4) = 1. 
For A(t) = the flow is steady and has hyperbolic stagnation points at: 



cos ( 



(35) 



hi = (0,0, -a) 
h 2 = (0,0, a), 



(36) 



where h\ has a one dimensional stable manifold and a two dimensional unstable manifold, 
and hi has a two dimensional stable manifold and a one dimensional dimensional manifold. 



In Branicki and Wiggins (2009) numerical algorithms are used to show that when A(t) 



increases from to the expression in (35), hi and h 2 "continue" to be hyperbolic (time- 



dependent) trajectories, denoted 71 and 72, respectively, and at each instant of time 71 has 
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Figure 16: a) The two dimensional unstable manifold of 71 is shown at t = 5.3 for U = a = 2 
(image courtesy of M. Branicki); b) the intersection of the unstable manifold with a vertical 
cross section of M2 obtained at y = —0.3 and t = 5.3 for r = 10; c) the intersection of the 
unstable manifold with a horizontal cross section of M2 obtained at z = and t = 5.3 for 
r = 10 



a one dimensional stable manifold and a two dimensional unstable manifolds, and 72 has 
a two dimensional stable manifold and a one dimensional stable manifold. The unstable 
manifold at t = 5.3 with U = 2 and a = 2 computed as in Branicki and Wiggins (2009) is 
displayed in Figure I61). 

We now illustrate the performance of Lagrangian descriptors for this three dimensional, 
aperiodically time dependent example. Figure [T6]o) and c) show respectively the intersec- 
tion of the unstable manifold with vertical (at y = —0.3) and horizontal (at z = 0) cross 
sections of M2 obtained under the same conditions. The agreement is excellent. 

The Lagrangian descriptor M2 provides many details on the Lagrangian structure of 
the flow as 17 1) confirms. It shows the results on the vertical section displayed in Figure 
[16)3) obtained with r = 10. Similar slices are displayed in Figure [17)3) for the descriptor 
M\. The output is similar to M2, but with a lower contrast which make more difficult the 
visualization. The structure is clean and sharp in contrast to that provided by methods 
discussed next. FTLE computations, similarly to what happened in the 2D case, distorts 
the image by introducing numerous structures that do not have Lagrangian interpretation. 
The output of forwards and backwards FTLE is displayed in figure 17 2) and d) confirms 



this extreme. Alternative averages as those proposed in section 2.3 provide the output 
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shown in figure 17 3) and f). These averages also present spurious structures that distort 
the output and blur the true Lagrangian information. 



4 Application of Lagrangian Descriptors to Velocity Fields 
Defined as Data Sets 

In this section we discuss the performance and application of Lagrangian descriptors to 
velocity fields defined as data sets. In particular, we will consider the finite time data set 
obtained from the numerical simulation of a wind driven, quasi-geostrophic (QG) 3-layer 



model in a rectangular basin geometry on a /3-plane (Rowley (1996)). The basic circulation 
pattern in the upper layer is an unsteady double-gyre structure, i.e. a counterclockwise 
gyre in the north and a clockwise gyre in the south separated by a jet. The velocity data 
set is obtained on a 1000 km x 2000 km rectangular domain with spatial and temporal 
resolutions of 12.5 km x 12.5 km and 2 hours (although the data is only saved every 24 
hours), respectively, and spans over 4000 days. The 4000 day interval of the data set is 
considered after the fluid is started from rest and allowed to spin up for 25000 days. More 
details of the specific parameters used and the numerical method used for solving the QG 



equations can be found in Rowley (1996) or Coulliette and Wiggins (2001) 



We have chosen this data set because it is aperiodic and results on the computation 
of distinguished hyperbolic trajectories (DHTs) and their stable and unstable manifolds 



have previously been reported in Mancho et al. (2004, 2006c); Madrid and Mancho (2009). 



In particular, we choose three hyperbolic trajectories that have previously been identified 



in a range of days between and 900 in Mancho et al. (2004, 2006c). Figure 18 shows 



one of these DHTs located in the northern gyre, which has been shown to possess the 
"distinguished property" in the interval between 120 days and 300 days (more details on 



this can be found in Madrid and Mancho (2009)). Beyond this time interval the trajectory 
still exists, but it no longer has the "distinguished" property. The loss of the distinguished 



property has been linked in Madrid and Mancho ( 2009 ) with the loss of hyperbolicity. We 



remark that in complex time varying geophysical flows one frequently observes trajectories 
whose finite time stability characteristics change from hyperbolic to elliptic, and vice versa. 
This is a particularly challenging phenomenon for the dynamical systems approach and has 



(2009); Branicki et al. (2011). 



been discussed in Branicki and Wiggins (2010); Mancho et al. (2006b); Madrid and Mancho 



Figures 19 1) and b) show the results of the direct computation of the stable manifold 



for two of the DHTs and the unstable manifold for the remaining DHT for days 290 and 



320 using the algorithms described in Mancho et al. ( 2004 ) . This figure also shows contour 
plots for a selection of Lagrangian descriptors computed for r = 150. 

We note that there is an essential difference between the figures showing the direct 
computation of DHTs and their stable and unstable manifolds and the contours of La- 
grangian descriptors. In order to compute manifolds of DHTs, we must first identify and 
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b) d) f) 

Figure 17: a) Evaluation of M2 on the vertical plane (x, y = —0.3, z) at t = 5.3 for r = 10; 
b) the structure provided by Mi also on the vertical plane (x, y = —0.3, z) at t = 5.3 for 
t = 10; c) the forward FTLE field on the vertical plane (x, y = —0.3, z) at t = 5.3 for 
t = 10; d) the backward FTLE on the vertical plane (x, y = —0.3, z) at t = 5.3 for r = 10; 
e) the forward average of the horizontal component of the velocity; f ) the backward average 
of the horizontal component of the velocity. 
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Figure 18: The time evolution of a DHT; a) the x-coordinate; b) the y-coordinate. Units 
are in thousands of kilometers. 



compute the DHTs. Lagrangian descriptors, on the other hand, compute the stable and 
unstable manifolds for all DHTs present during a chosen time interval, i.e., there is no 
need to identify the DHTs a priori. 

Typically the Lagrangian structure revealed by different Lagrangian tools contains more 
details for longer integration times r. However, it is also of much interest in applications 
to obtain accurate Lagrangian information for small r. We will address how Lagrangian 
descriptors perform in this context, as well as make comparisons with the performance 
of FTLEs, by studying the convergence rate of different Lagrangian descriptors with r 
towards the singular lines that reveal the Lagrangian structure of the flow. Our attention 
is not on the manifold structure for a large spatial region, but on the stable and unstable 
subspaces tangent to the stable and unstable manifolds of a particular DHT. We remark 



that these subspaces are related to the also called Lyapunov vectors Legras and Vautard 



(1996) that indicated dispersive directions in the fluid flow, backwards and forwards in 



time. 

We place our study at day 130, far from the transition of the DHT. Figure [20| shows the 
stream function on this date with the position of the DHT indicated at (70.631, 1357.660) 
km. We use Lagrangian descriptors to approximate the stable and unstable subspaces of 
the DHT as follows. We consider a "small" circle centred at the DHT. In particular, we 
consider the circle centred at the DHT parametrized as (r cos(a), r sin(a)) and of radius 
r = 1.5 km. We then compute the Lagrangian descriptor for trajectories starting in this 
circle for different r, and consider the singular contours of the Lagrangian descriptor near 
the DHT. These are approximations to the stable and unstable subspaces of the DHT. 



Figure 21 summarizes the convergence of different Lagrangian descriptors towards the 
stable and unstable subspaces of the DHT at day 130. The angle positions of the stable 
and unstable subspaces are marked with vertical lines. Panel 21 1) shows Mi, M3 and M5 
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e) * f) * 

Figure 19: Pieces of stable manifold (blue) for two DHTs and and unstable manifolds 
(red) of the remaining DHT shown on a) day 290; b) day 320. Contour plots for a selection 
of Lagrangian descriptors evaluated over the quasigeostrophic flow with r = 150. c) M3 at 
day 290; d) Mi at day 320; e) M 5 at day 290; f) M 2 at day 320. The units on the axis are 
in thousands of kilometers. 
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X 



Figure 20: Stream function of the quasigeostrophic model at day 130. The position of a 
DHT at this day is marked with a circle. The units on the axis are in kilometers. 



versus a for r = 45, which for this specific vector field is not a large value. At this stage 
singular features are only developed by M3. The Lagrangian descriptor Mi shows a slower 
convergence towards the subspace aligned along a = 2.032 and a = 2.032 + ir/2. The 
features displayed by M5 are not yet sharp enough at this stage. Panel [2T]b) confirms a 
finer structure for r = 65 for the Lagrangian descriptor M3. The Lagrangian descriptors 
M4 and M5 perform worse at this stage, while the output of Mi and M2 is similar to that 



of M3. Panel 21 3) shows M2 and M4 at r = 95, when convergence has been reached for all 
of the Lagrangian descriptors considered. 

These results indicate that the Lagrangian descriptor M3 performs the best in the sense 
that it indicates the subspaces for smaller r. It is followed by Mi, M2, and M5, respectively. 
The Lagrangian descriptor M4, also develop sharp features along the same structures, but 
it require longer time intervals of integration for the situations considered here. 

Figure 22 contrasts the ability of FTLE to locate the stable and unstable subspaces of 
the DHT with that of the Lagrangian descriptors. For comparison purposes, we choose M3. 
For reference, the angle positions of the stable and unstable subspaces are marked with 
vertical lines. Panel [22^i) shows the forward and backward FTLE and M3 versus a for r = 
45. The Lagrangian descriptor M3 is represented by a thick black line. Its sharp features 
mark accurately the position of the stable and unstable subspaces. The forward FTLE, 
represented by the thin black line, at this stage presents pronounced features, although not 
yet sharp. Deviations from the correct position of the stable subspace are significant, and 
the performance of forward FTLE is clearly lower than that of M3. The backward FTLE, 
represented by the thick gray line, accurately detects the unstable subspace with a similar 
performance to that of M3. Panel [22|j) presents results for r = 65. However, deviations 
of the forward FTLE field (thin black line) with respect to the stable subspace are still 
noticeable in this panel. Panel [22p) confirms that at r = 95 both FTLE and Lagrangian 
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b) 




Figure 21: a) Lagrangian descriptors M\ (black), M3 (red) and 100 x M5 (blue) versus a 
for r = 45; b) descriptors 0.2 x M3 (black), 10 x M4 (red) and 10 x M5 (blue) versus a for 
r = 65; c) M2 (blue) and M4 (red) versus a for r = 95 
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Figure 22: a) Forwards FTLE (black), backwards FTLE (red) and M3 (blue) versus a 
for r = 45; b) forwards FTLE (black), backwards FTLE (red) and M3 (blue) versus a 
for r = 65; c) forwards FTLE (black), backwards FTLE (red) and M3 (blue) versus a for 
r = 95 



descriptors perform well at locating the angles of the stable and unstable subspaces. 



5 Computational performance 



Computationally the evaluation of Lagrangian descriptors offer numerous advantages in 
comparison with the computation of FTLE, as we discuss next. Calculation of FTLE 
requires the evaluation of the gradient of the flow map, which is numerically realized by 
expression (49) in the Appendix. As indicated there for obtaining the a field at a grid 
point x* 7 - in a 2D flow it is requested to integrate the dynamics of four neighbour points 

in order to keep valid 
Thus in order to be accurate, the 



hi 

-lj> x i+l 



,x 



1,3- 



_ 1; x* - +1 which must be a small distance e from x.* j 



expression (49) in the range of integration (t*,t* + r) 



numerical approach (49) must keep a balance between the distance of grid points e and 



the integration time r. These parameters must be appropriately tuned. This tuning is not 
necessary in the computation of the Lagrangian descriptors, thus making their computation 
more direct. On the other hand if the grid size in which the a field is evaluated is bigger 
than e, then for obtaining its value in a point x*^ expression (49) indicates that four 
integrations are necessary, while Lagrangian descriptors require only one. This makes 
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Lagrangian descriptors computationally much more efficient, specially in geophysical flows 
given as finite time datasets, because in them integrations require interpolations which are 
computationally expensive. Some approaches compute FTLE taking e as the grid size in 
which the a field is evaluated. In this case integrations of neighbours are stored to prevent 
multiple performaces of the same trajectory, and then only one integration per point x* • 
is requested. However this approach has also several computational disadvantages with 
respect to the use of Lagrangian descriptors. One is that it is harder to program because 
the evaluation of the FTLE at each point requires access to stored information, and the 
second is that the grid size is doomed to be small enough to keep expression ( 49 ) accurate 
in the range (t*,t* + r). This size can be much smaller that the grid size required to have 
a smooth representation for the lagrangian descriptor, and in this case the computation of 
FTLE is again more expensive. 

Another difference between Lagrangian descriptors and FTLE is that the former provide 
the stable and unstable manifolds simultaneously in the same output, while the latter 
require post-processing of the results into one picture. In the post-processing of the output 



typically the results are 'filtered' by using a threshold in the field values (see Branicki and 



Wiggins (2010)) and in this processes spurious structures may be removed. In the case 



of Lagrangian descriptors an accurate image of the manifold features is obtained without 
the need of any 'filtering' process, de la Camara et al. (2012) have demonstrated that in 
applications having both the stable and unstable manifolds simultaneously in one picture 
is an interesting feature that is essential for the successful analysis of transport events. 



6 Conclusions and Outlook 

In this paper we have proposed a general method for revealing geometrical structures 
in phase space that is valid for aperiodically time-dependent dynamical systems. The 
central quantities in our method are referred to as Lagrangian descriptors and are based 
on finite time integrals of positive, bounded intrinsic geometrical or physical properties 
of the dynamical system along trajectories. In particular, these properties are integrated 
forward and backward in time on some interval (t — r, t + r). In this way, the Lagrangian 
descriptors have the capability of revealing both the stable and unstable manifolds in the 
same calculation. We have shown that the convergence to geometrical structures of the 
dynamical system requires the use of a long enough r. For r values below a threshold, which 
depends on each particular dynamical system, the pattern provided by the Lagrangian 
descriptors is far from the geometrical structures. However beyond that threshold, for 
increasing r, the output of the Lagrangian descriptors becomes more and more complex, 
revealing more details of the geometrical structures of the dynamical system. We have 
given an analytical justification for these observations for the benchmark case of a linear 
vector field having a hyperbolic saddle point. 

In all examples that we have considered we have carried out simultaneous comparisons 
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of FTLEs and finite time averages of components of the vector field along trajectories 
with Lagrangian descriptors. We have shown that Lagrangian descriptors provide superior 
performance in all cases, in the sense that they more accurately reveal the geometrical 
structures and they require a shorter time to converge to the geometrical structures. For 
example, we have shown that FTLE give spurious structures for the integrable Duffing 
equation and fail completely to provide any structure for the linear saddle. Finite time 
averages of the horizontal component of the vector field exhibit the same deficiencies as 
FTLEs for these particular two examples. We have also shown that Lagrangian descriptors 
accurately provide the geometrical structures for a three dimensional flow. We have com- 
pared these results to results obtained from FTLE calculations and finite time averages 
of a component of the vector field along trajectories and we have seen that Lagrangian 
descriptors provide a more accurate representation of the geometrical features of the flow 
than these techniques provide. We have examined the issue of the convergence time to the 
stable and unstable subspaces of a DHT for a velocity field given as a data set and have 
shown that Lagrangian descriptors converge more rapidly to these subspaces than FTLE. 
This is promising from the applications point of view, in which reductions for the fore- 
cast time of the velocity field is desirable. Lagrangian descriptors can give more accurate 
predictions with less information in time. Finally we have shown that computationally 
Lagrangian descriptors provide several advantages over FTLE. In general, their use leads 
to a reduction in CPU requirements by a factor of 4, their implementation in code is sim- 
pler, and they require no "parameter tuning". Moreover, Lagrangian descriptors provide 
both the stable and unstable manifolds in one output. In order to achieve this with FTLE 
post-processing of the stable and unstable is required in order to combine the two outputs 
into one picture at the same time. 
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A Numerical computation of Lagrangian descriptors 



In this appendix we describe the numerical computation of general Lagrangian descriptors: 

M(x*,i*) v>T = f + F(x(t))dt. (37) 

Jt*-T 

where J r (x(i)) is a positive bounded scalar representing a geometrical or physical property 
that is evaluated along a trajectory. This scalar properties introduced in this paper depend 
on vector quantities such as velocity, acceleration, the time derivative of acceleration or 
combinations of these quantities (e.g. as in the example using curvature). For example, for 
the Lagrangian descriptor M±, Eq. ([2]), is evaluated by integrating the function ||v(x(t), t)|| 
along a trajectory, x(t), x(t*) = x* from from t* — r to t* + r. 

Geometrically, the computation of Mi is equivalent to evaluating the area A below the 
graph I |v(x(t), t)| | in the specified time interval. The area A is obtained from the following 
one dimensional dynamical system: 

dY 

— = (38) 

For the initial condition Y{t*) = 0, the area A is given by the value of Y at t* + r minus 
the value of Y at t* - r , i.e., Y(t* + r) - Y(t* - r) = A. 



The integration of the system (38) is performed with a 5th order variable time step 



Runge-Kutta method, in particular with the subroutine rkqs described in Press et al 



(1999). A peculiarity of this differential equation is that it depends on t both explicitly 
and implicitly through the trajectory of ([I]). A Runge-Kutta step from to to t\ applied to 
Eq. (38) requires evaluating ||v|| along trajectories of ([!]) at intermediate steps to + At. To 
this end the argument x(t) that must be passed to ||v|| at time to + At must be obtained by 
evolving the trajectory of ([!]) from (to, x(to)) to (to + At, x(to + At)). We remark that this 



adaptive method can be faster than the method proposed inlMadrid and Mancho (2009) 



which is based on fixed time step integrations. Moreover, this method is quite versatile, 



since from one descriptor to another it is only the right hand side in Eq. ( 38 ) which needs 
to be modified. 

In order to evaluate the entire family of descriptors presented in this paper three dif- 
ferent types of vectors must be obtained: velocity, acceleration and the time derivative of 
acceleration. Once these are obtained individually, their combinations relevant to the par- 
ticular type of descriptor of interest are straightforwardly computed. Velocity is a vector 
which is directly provided if Eqs. are given; the others are obtained as we explain next. 
The general expression for a(x, t) is as follows: 



dv dv dv dv 

dt dt dx x dy 



a ( x > *) = = a7 + nZ v x + ^7 v y (39) 
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which is easily computed if there is an exact expression for the velocity field v. If the 



velocity field is given as a finite time data set, Eq. (39 ) is evaluated by providing a discrete 
version of the equation: 



as follows: 



dv d 2 x 

~dt ~ ~d¥ 



a(x(t), t) = W + h),t + h)-vW-h),t-h) + 0(h2) 



The time derivative of the acceleration is given by: 



(40) 



(41) 



da(x, t) 
dt 



d dv 

dt dt 

dxdt 

d 2 v 



d 2 



cPv 

dt 2 dtdx 



v dv dv^ 



5 2 A 



dv dv„ 



dx dt dtdy dy dt 



+ 



+ 



tfv 

dx 2 

Q2, 



V x + 



dv dv q 



+ 



d 2 v 



dx dx dxdy 



v v + 



dv dv y 
dy dx 
dv dv 



+ 



_ "v dv dv x d 2 v 
dydt dydx x dx dy dy 2 y dy dy 



(42) 



If the velocity is supplied as a finite time data set, Eq. (42) is evaluated by providing a 
discrete version to the equation: 



da 

~dt 



d^v 
dt 2 ' 



(43) 



which reads as follows: 



da(x(t), t) v(x(i + h),t + h) - 2v(x(t), t) + v(x(t -h),t- h) 



dt 



2h 



+ 0(h 2 ). (44) 



For the results reported in this article Eqs. (41) and (44) are implemented where needed 
using h = 10 -4 . 



B Computation of FTLE's 



In this appendix we give a brief discussion of the computation of FTLEs, following Lekien 
et al. (2007). The FTLE field is denoted by cr^,(x*), and is computed from the maximum 



eigenvalue A max of the matrix A: 

4(x*)= 1 ln V / WA)> (45) 



T 
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where the matrix A is defined from the flow map, as we now describe. Consider a point 
x* at time t* , and let it be advected by the flow map, (f>, from t* to t* + t: 



x* -+ 4: +t (x*). 

The gradient of the flow map defines the matrix N: 

"(X*) 



N 



dx* 



(46) 



(47) 



Using the convention that N T denotes the transpose of N, A is the symmetric matrix: 

A = N T N (48) 

Note that in two dimensions the gradient of the flow map at a grid point x* ■ = (xij,yij) 
can be computed central differences as follows (see Shadden (2005)): 



(x*) 



(fx* 



Xi + lj(t*+T)-Xj_lj(t*+T) XjJ + l(t*+T)—Xjj-l(t*+T) 

Xi+i,j(t*)-Xi-ij(t*) Vi,i+i(t*)— l/t,J-l(**) 

yi+l,j(t*+T~)-yi-i,j(t*+T) yij + i(t*+r)-j/i,j_i(t*+r) 



(49) 



where (^-(i* + r), y fi (f + r)) = 4>\+ T {xUt*)). 



C Regularity of Lagrangian descriptors 

In this appendix we discuss some issues associated with the regularity of Lagrangian de- 
scriptors, especially when dealing with velocity fields defined as data sets, as discussed in 
Section |U 

In this case, since the data is discrete in space and time, in order to have a continuous 
representation for the dynamical system, 



dx 



v(x,t), x <E 



(50) 



the velocity field v must be interpolated in space and time in order to compute particle 



trajectories. In Mancho et al. (2006a) it was shown that bicubic interpolation in space 



and 3rd order Lagrange polynomials in time is sufficiently accurate for our purposes. This 
interpolation is C 1 continuous in space and C° continuous in time. 

Lagrangian descriptors are sensitive to the quality of the interpolation of the data field, 
as are FTLE. Regarding Lagrangian descriptors, we focus first on those involving the time 



derivative of the acceleration. Following Eq. (42), the time derivative of acceleration along 



trajectories involves second order spatial derivatives of a velocity field (50) which are not 



so regular. Figure 23 1) illustrates the evaluation of the time derivative of acceleration 
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Figure 23: a) Evaluation of the x component of the vector da/dt along a trajectory from 
t — t to t + r; b) evolution of M4 in the same time interval; c) evolution of the time 
derivative of M4. 




Figure 24: Evaluation of the x component of the vector a along a trajectory from t — r 
to t + r; b) evolution of M2 in the same time interval; c) evolution of the time derivative 
of M 2 . 



along a trajectory for the system (50) with initial condition t = 130 days, x = 410 km, y = 
1100 km. The integration is performed backwards from t to t — r and forwards from t 
to t + r, with r = 20. In particular, this figure shows the evolution of the horizontal 
component of this vector. 



Despite the lack of regularity displayed in Figure 23 1) , the actual descriptor is more regular 
since it is based on integrals over these functions. Figure [23|j) shows the evolution versus 
t of M4. Figure 23:) displays the time derivative of M4, confirming that its regularity 



is comes from that of da/dt. Expression (39) shows that the evaluation of acceleration 



along trajectories involves first order spatial derivatives of the velocity field (50). Figure 



24 1) shows the evolution of the horizontal component of acceleration, which confirms that 
the regularity of the interpolator of the underlying velocity field is observable. The actual 
regularity of M2 is larger than that of acceleration as it is based on integrals over these 
quantities. Figure [24)j) shows the evolution versus r of M2 as a more regular curve, and 
in c) the time derivative of M2, confirms that its regularity comes from that of a. 
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Figure 25: a) Evolution of d y v x (x(t)) versus t for backwards intergation; b) evolution of 
the N\ t 2 component versus t for backwards integration; c) evolution of the time derivative 
of the backwards FTLE. 



The performance of FTLE is similar to the performance of Lagrangian descriptors in- 
volving the evaluation of acceleration, as they require the evaluation of the gradient of the 
flow map. The gradient of the flow map is based on integrations on linear equations which 
are obtained by linearizing the velocity field (50) along a trajectory. This involves evalu- 
ating first order spatial derivatives, as in the case of acceleration. Figure 25 1) shows the 
partial derivative of v x with respect to y along the same trajectory as above, for backwards 
integration. The lack of regularity is noticed, however as the FTLE computation is based 
on integrations along this linearized field, it gains regularity and for instance in Figure 
^5p) we show the evolution of one of the components of the matrix N used to compute 



the FTLE. Again in Figure 25;) illustrates how the regularity of the backwards FTLE is 
governed by the regularity of linearized velocity fields. Finally we note that descriptors 
based on the integration of the velocity are the most regular. They are advantageous in 
this regard not only to those LD involving acceleration or its time derivative, but also to 
FTLE. 

The discussed reasons suggest that for analyzing velocities fields given as data sets, 
the choice of Lagrangian descriptors involving \da/dt\ may be less appropriate than those 
involving v or a since they require interpolators with a higher order of regularity than is 
required by the latter quantities. On the other hand, for data sets interpolated with not 
too regular interpolators Lagrangian descriptors based on v are the better choice, even 
with respect to FTLE. 



51 



